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Abstract. We demonstrate the application of transition state theory to wave packet 
dynamics in metastable Schrodinger systems which are approached by means of a 
variational ansatz for the wave function and whose dynamics is described within the 
framework of a time-dependent variational principle. The application of classical 
transition state theory, which requires knowledge of a classical Hamilton function, 
is made possible by mapping the variational parameters to classical phase space 
coordinates and constructing an appropriate Hamiltonian in action variables. This 
mapping, which is performed by a normal form expansion of the equations of motion 
and an additional adaptation to the energy functional, as well as the requirements to 
the variational ansatz are discussed in detail. The applicability of the procedure is 
demonstrated for a cubic model potential for which we calculate thermal decay rates 
of a frozen Gaussian wave function. The decay rate obtained with a narrow trial 
wave function agrees perfectly with the results using the classical normal form of the 
corresponding point particle. The results with a broader trial wave function go even 
beyond the classical approach, i.e., they agree with those using the quantum normal 
form. The method presented here will be applied to Bosc-Einstein condensates in the 
following paper [A. Junginger, M. Dorwarth, J. Main, and G. Wunner, following paper, 
submitted to J. Phys. A]. 
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1. Introduction 

Transition state theory (TST) is one of the most important tools for analysing chemical 
reactions qualitatively as well as quantitatively and has also applications remote from 
its origin, e.g. in atomic and semiconductor physics [l|[2], the study of clusters [S], and 
cosmology [2]. 

Fundamental to all these calculations is a classical Hamilton function H{q,p) given 
in phase space coordinates q, p which is either an inherent part of the underlying model 
or, e.g. for chemical reactions, given by an ab initio calculated potential energy surface 
V{q) 

In this paper and the following |9j, we apply TST to wave packets, which are widely 
used for variational approaches to quantum mechanical systems, such as Bose-Einstein 



condensates 10 ,11 . Here, the Schrodinger or Gross-Pitaevskii equation is solved within 
the Hilbert subspace of the variational parameters, and to describe the dynamics of the 
wave function, one makes use of a time-dependent variational principle (TDVP). 

In order to apply TST to a system with a known Hamilton function, one transforms 



the latter to its normal form using canonical transformations 12 . In contrast, for 
the application of TST to wave packets, where such a Hamiltonian is not known, we 
demonstrate a method to locally map the variational parameters to action variables 
in classical phase space using non-canonical transformations. This step is performed 



by a normal form expansion 13 of the equations of motion derived from the TDVP 
and the corresponding transformation of the energy functional. After having expanded 
the equations of motion into normal form, they are integrated to a common Hamilton 
function and an additional transformation guarantees its equivalence to the energy 
functional, which at the end serves as the classical Hamiltonian. 

To demonstrate the procedure, we apply this method to calculating the decay rate 
of a frozen Gaussian wave function located at the metastable ground state of a cubic 
model potential and compare it with the ones resulting from ordinary classical and 
quantum normal forms |12j of the corresponding point particle. For the application to 
Bose-Einstein condensates with coupled Gaussian wave functions, we refer the reader to 
reference [o]. 

Our paper is organized as follows: First, we give a brief review of the classical and 
quantum normal form expansions, which can be applied to a given Hamiltonian H{q,p). 
Then, we introduce the variational approach, demonstrate the general concept of 
mapping the variational parameters to classical phase space and illustrate the calculation 
of the respective thermal decay rates. At the end, we apply the three methods to the 
model potential and compare their results. 
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2. Theory 

2.1. Classical and quantum normal form 

We begin by summarizing the classical and quantum normal form transformation of a 
given Hamilton function 

Hiq,p)=Tip) + V{q) (1) 

in the vicinity of an equilibrium point Qq. Since this is a well known theory, we will 



only give an overview of the main aspects and refer the reader to Waalkens et al 12 
and references therein for details. 

In order to obtain the normal form Hamiltonian, one expands the Hamilton function 
g into its power series around the fixed point Qq, 

oo 

H = Y,Hn. (2) 

n=0 

Here, Hn summarizes monomials which are homogeneous of degree n, and for simplicity, 
we omit the dependency of the phase space variables q,p from now on. 

After having shifted the fixed point to the origin and having applied a symplectic 
transformation which "simplifies" the quadratic part H2 in a way that diagonalizes the 
corresponding Hamiltonian equations of motions, the Hamilton function is in classical 
normal form with respect to the quadratic part. The aim is to transform the whole 



Hamilton function H into normal form, which is defined by the condition 12 

adjy.i/ = {H2, H} = (3) 

and where the adjoint operator ad/^j-f^ = {H2, H} is defined by the Poisson bracket. 

To transform the higher order terms n > 3, we assume the Hamiltonian to be in 
normal form up to order n — 1 which we denote by the superscript and apply 

successive symplectic transformations given by the generator Wn- This generator Wn, 
which is a solution of the homological equation 

^(n) = ^{n-l)_|^(2)^^^|^ (4) 

then transforms the Hamilton function according to 

00 _ 

fc=0 

These steps are performed successively for each order n until the Hamiltonian is in 
classical normal form i^cNF up to the desired order. 

In analogy to the classical version a quantum normal form can be calculated by 
replacing the adjoint operator adw by the Moyal adjoint action 

2„, . 1/ , 



12 



14 



%X)-(%X)]]h. (6) 



MadwH = - Wsin . 

h \2 

Since the procedure of getting the quantum normal form is very similar to that of its 
classical analogue, we skip the details here and again refer the reader to references 
T2lfT4 . 
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Altogether, for a one degree of freedom Hamilton function of the form 

H = ^^ + V{q\ (7) 

with stationary point go we will use it for comparisons with the formalism for wave 
packets in section [3| the quantum normal form Hamiltonian Hq^y can explicitly be 
written down in terms of the action variable J, and in order 0{J^) its symbol reads 



12 





1152A 



Here A = \/—V"{Q) is the eigenvalue of the linearised Hamiltonian equations of motion 
and all derivatives are evaluated at q = Qq. The classical normal form Hqtssf of the 
Hamiltonian ([T]) can easily be obtained from equation ([s]) by setting /i — > 0. 

2.2. Variational approach 

We now focus on the system from a different point of view and describe it quantum 
mechanically via its wave function ip{r,t). Its physical properties are given by the 
Hamilton operator H and the time evolution of the wave function is determined by the 
time- dependent Schrodinger equation 

Hijir,t) = idtij{r,t), (9) 

where we set h = 1 from now on. 

In the following, we will approach the system in the framework of a variational 
ansatz 

t/j{r,t)=^{r,z{t)) (10) 

for the wave function in the form of a wave packet. Here z{t) = {zi, . . . , ZdY" is a set 
of d complex and time-dependent variational parameters and the time evolution of the 
wave function is fully determined by that of the parameters. 

A usual way of treating such a variational ansatz is by means of a variational 
principle 15 applied to the time-dependent Schrodinger equation (|9]). However, we 



will also discuss it from the perspective of field theory with the Hamiltonian density 
T-i. Although the latter approach will give exactly the same equations of motion for 
the variational parameters and the same energy functional, it allows for a better insight 
into the "Hamiltonian relation" between those two, and is the basic foundation of the 
variational approach and the subsequently performed normal form expansion. 
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2.2.1. The time- dependent variational principle An approximate solution of the time- 
dependent Schrodinger equation (|9]) within the parameter subspace of the wave function 



(10) is given by the McLachlan variational principle 15 
/ = ||i0 — H''p\\ = min. 



:ii: 



where the quantity / is minimized with respect to and (p = ip is set afterwards. For 



the parametrized wave function (10) the application of the TDVP results in a set of 



first-order differential equations 16 
Kz = -ih. 



(12) 



which determines the time evolution of the variational parameters, and where the matrix 
K and the vector h are defined by (m, n = 1, . . . ,d) 

d-ip \ * dip 



dV 



dV 



dip 



Hip. 



(13) 
(14) 



The energy of the system, which is given by the expectation value 

E{z) = Jd'r {ip{r,z{t))r Hip{r,z{t)), 
is also a function of the variational parameters. 



(15) 



2.2.2. 



Approach from field theory Within the McLachlan variational principle. 



equation (11), one recognizes no obvious "Hamiltonian relation" between the energy 



functional (15) and the equations of motion (12). By contrast, such a relation becomes 



obvious, if we focus on the ansatz from field theory with the Hamiltonian density 
"H = ip*Hip. The energy functional is then given by 



Eiz) 



(16) 



which is precisely the result of equation (15). 



Within field theory, the dynamics of the wave function ip is given by (with the field 
momentum tt = iip*) 

m 



ip 



idip* 



which, for the variational ansatz ( 10 ) of the wave function, leads to 

dip 



dip 
dz„. 



dz* 



(17) 



(18) 



and, after integration over d^r, is exactly the result of equations (12)-(14). This step 



only needs a few lines of elementary calculations, including ip = Q^Zn, d^* 



dip- 



a multiplication with the inverse of the derivative [dz^/ dip*]~^ = dip* /dz^ and the 
two assumptions that dip*/dz'^ = (dip/dzm)* and that the Zm be complex, i.e. the 
derivatives with respect to z^ and z^ are independent. Moreover, we note that for ip* 
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the corresponding equation of motion ip* = —dT-L/idi/j yields the complex conjugate of 



equation (18). 



Vice versa, from equation ( 18 ) we see that the equations of motion derived from the 



TDVP do not only hold in the "integrated version", equations (12)-(14), but in these 
equations, we may set the integrands equal to one another, and those directly follow 
from the Hamilton operator H and its corresponding Hamiltonian density H. 

From this point of view, a "Hamiltonian relation" between the energy functional 



(15) and the equations of motion (12) is obvious, and in addition, we obtain the 
conditions for the variational parameters to be complex and the variational ansatz for 
the wave function to fulfil dip* /dz^ = {dip / dz^)* ■ The latter requirement directly 
follows from the Cauchy-Riemann differential equations, i.e. the variational ansatz for 
the wave function has to be complex differentiable in the variational parameters in 
order to obey this relation. 

2.3. Mapping to phase space 

Within the variational approach, the dynamics of the wave function and, with it, any 
"reaction" in the system, is given by the vector z of the complex and time-dependent 
variational parameters. TST may then be applied to the wave packet if one is able to 
divide the space of the variational parameters - after a suitable definition of the latter - 
into "reactants" and "products" , and if the point with the least energy on the dividing 
surface which forms the "activated complex" of the system is a stationary point of 



saddle-centre-. . .-centre type of the energy functional (15). The corresponding reaction 
rate is then qualitatively given by the flux over this saddle in the complex z-space. 
However, in order to calculate the flux quantitatively a corresponding Hamiltonian in 
phase space is required. 

The definition of the dividing surface and the calculation of the flux were given in 



reference 12 , for the case that the corresponding Hamiltonian of the system is known. 



This is, however, not the case within the variational approach which is based on the 



energy functional (15) and the equations of motion (12). A globally defined Hamilton 



function which describes the system is not known, in particular the energy functional 



(15) cannot serve as a Hamiltonian because of the fact that the equations of motion (12) 
cannot be derived from the former via the canonical equations, and canonical variables 
are not even defined therein. 

Nevertheless, as will be shown below it is possible to construct a local Hamiltonian 
by an appropriate change of variables from the variational parameters z to local 
phase space variables p, q which equivalently describes the dynamics of the variational 
parameters and at the same time reproduces the energy of the system, given by equation 



(15). Therefore, this constructed Hamiltonian can be applied to calculating the flux. 

In the following subsections, we provide the general concept of how to construct 
this equivalent Hamilton function. It consists of the subsequent steps: 



(i) Taylor expand the equations of motion (12) in the vicinity of a fixed point Zq in 
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terms of the variational parameters and diagonalize the latter with respect to their 
linear part. 

Perform a normal form transformation of the Taylor expanded equations of motion 
and introduce canonical variables. 



(iii) Analogously to step (i), Taylor expand the energy functional (15) and apply the 
change of coordinates which corresponds to the transformation from step (ii). 

(iv) Integrate the equations of motion according to Hamilton's equations and adapt the 
resulting Hamiltonian H to the transformed energy functional to obtain the final 
Hamilton function H. 



2.3.1. Transformation of the vector field We assume the differential equations (12) 
to exhibit a fixed point Zq. To obtain the local classical Hamilton function, we first 
Taylor expand these equations of motion in the vicinity of the fixed point up to the 
order nmax and split the expansion into its real and imaginary part. This results in the 
2(i-dimensional real vector field 

X = a{x) = CLn{x) (19) 
ra=l 

where x = {Re{zi—zoi), Im(zi — zqi), • • • , ^^izd — zod), l^{zd — ZQd)) denotes the deviation 
of the variational parameters from the fixed point and a„(a;) summarizes all terms of 
the expansion which are homogeneous of degree n. 

In order to "simplify" these differential equations we diagonalize them with respect 
to their linear part ai{x) = Ai-x in the next step. Since this is trivial, we will not go into 
it in more detail. Note, however, that, even if it does not affect the result of the linearised 
and diagonalized equations of motion, this procedure itself is not unique, because one 
may work with arbitrary complex multiples of the eigenvectors of which one makes use 
of in this step. In contrast to the linearised differential equations, the choice of these 
multiples does affect higher order terms and the corresponding transformation of the 



energy functional which is explained in section 2.3.2 We will use this freedom, among 
others, to adapt the "integrated Hamiltonian" to the transformed energy functional in 
a last step. 

To further "simplify" the higher order terms of the remaining set of differential 
equations, we perform a near-identity transformation x ^ y given by (cf. reference 
[13]) 

X = cj)^{y), X = 0^=0(2/) = y (20) 
where the parameter e is chosen in such a way that we obtain the identity-transformation 
for e = and x = 4>^{y) solves the differential equation 

t-oi^) (21) 

with g{x) being the generating function. Concerning the normal form transformations 
which will be applied in the following, we keep close to the notation introduced in 



reference 13 and also use the notion "normal form" as defined therein. 
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Under the change of variables (20), given by a specific generating function g{x) the 

(22) 



vector field (19) transforms as [13 



oo ^ 



y 



k=0 



x=y 



with the Lie-operator Cg defined by 



Cga{x) 



g{x) ■ V 



a[x] 



a(x) ■ V 



9[x)- 



(23) 



By simply expanding the sum in equation (22 ) and sorting by the polynomial order 



n it can, with b{y) = b^iy), easily be shown that a generating function 5'„(a;) which 
is homogeneous of degree n transforms the corresponding term an{y) according to the 
homological equation 

Kiy) = arXy) + Cg^a^{y). (24) 

Terms of CL{y) of lower degree than n remain unchanged. Therefore, by choosing a 
specific degree n of the generator g^, we can normal transform the vector field ci{y) 
order by order. 

Introducing the expansions 



[y) 



\m\=n 
\m\=n 



(26) 
(27) 
and \m\ = J2i 

with the integer vector m, and inserting these into equation (24) the transformation can 



9niy) = Yl ^mV"^^ 

\m\=n 

where we made use of the multi-index notations x'' 



Xi ^2 J'S 



be directly applied to the single coefficients which componentwise transform according 
to 

f3mi = Oimi + (Am - \i)-fmi ■ (28) 

Here, the Aj are the eigenvalues of the linearised equations of motion. In case of 
Am— Aj 7^ the respective monomial ami can be eliminated {f3mi = 0) by an appropriate 
choice of •jmi whereas in the "resonant" case 

Am - Ai = (29) 

it remains unchanged {Pmi = otmi) and cannot be eliminated. Hence, the eigenvalues 
of the linearised equations of motion determine the whole structure of the normal form 
expansion. 



Since the equations of motion derived from the TDVP, equations (12)-(14), are 



invariant under time reversal t — t- — t (together with {ip,z) — )■ {ip*,z*)), the eigenvalues 
of the linearised equations of motion in general exhibit the structure 

Xi = ±6i, . . . , ±6j, ±iwj+i, . . . , ±iu}d (30) 
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with 5i,uji G R, i.e. they occur pairwise with different sign. Assuming rational 
independence of the eigenvalues 5j,iwj {i = 1,2, ... ,d) this structure directly implies 



that equation (29) is fulfilled for 



m2i-i=m2i + l, {i = 1,2, . . . ,d), (31) 
m2j-i = m2j, (j 7^ i). (32) 

Consequently, the form of the normal form transformed equations of motion is fixed 



because only monomials satisfying equations (31 ) and (32) will remain. We now assume 



these equations of motion to be sorted in a way that the linear terms occur blockwise 



according to the eigenvalue structure (30) where each block X2i~i,i2i ij- = 1, • • • , d) 
contains the same eigenvalue differing only in its sign. Then they take the form 

^(^2,-lX2,)'"^^ (33) 



/-'m(2j-l)J'2j-l -^21 



X2i-1 

X2i =5^/3m(2.) a;^i7^"'^rn(^2,-iX2,)™^^ (34) 
We now systematically introduce "canonical" coordinates = X2i-i and momenta 



Pi = X2i {i = 1, . . . , d) with the help of which we can rewrite equations (33) and (34) as 



Pi 



Pm(2i^l)qi p. 



2j 



(35) 
(36) 



2.3.2. Transformation of the energy functional To obtain the local energy functional in 
phase space coordinates q,p, we need to apply the change of coordinates corresponding 
to the normal form transformation of the equations of motion also to the expanded 
energy functional 

^maxH" 1 

E= ^n{x). (37) 



n=0 



Since the generating function is, for each order, known from the normal form 
transformation of the vector field ci{x), the explicit change of coordinates corresponding 



to the transformation (24) can easily be obtained. In dependence of the generating 
function it reads |13l 



y 



oo ^ 

y-v 



X, 



where the "right-hand multiplication operator" Vg is defined by its action 
on differentiable functions f{x). 



(38) 



(39) 



TST for wave packet dynamics. I. Thermal decay in metastable Schrddinger systemslQ 



Inserting the transformation corresponding to each step into the expanded energy 



functional (37) yields the form 



(40) 



The fact that the energy functional can be expressed in products qjPj, is here due to 
the "Hamilton-like" relation between the equations of motion derived from the TDVP 



and the energy functional, explained in section 2.2.2 



2.3.3. Construction of the Hamilton function H The exponents of the qi,Pi in the z-th 



component of the transformed equations of motion (35)-(36) exactly differ by one and 



do not differ for qj.Pj (j ^ i) which is, both, due to the eigenvalue structure (30) and 
leads to the fact that they can easily be integrated to a common Hamilton function 



H 



m2i 



■{qiPi 



m2. 



m2j 



according to Hamilton's equations 

_dH _ dH 

This is possible, if the remaining coefficients (3^, 
satisfy the conditions of integrability 

/5m(2i-l) 



for all 



/3m(2i) 
: 1,... 



d 



- ~Pm{2i)) 
Pm{2i'-1) 

m2i' 

Pm{2i') 

m2i'-i 



(41) 



(42) 

after the normal form expansion, 

(43) 
(44) 

(45) 



Note that equation (43) is automatically fulfilled 



if the diagonalisation of the equations of motion is performed with, in each case, 
two corresponding complex conjugate eigenvectors which are normalized identically. 



However, the equations (44)-(45) are, in general, not automatically satisfied after the 



transformations performed above. 

In a system with one degree of freedom {d 
and the formal integration to the Hamiltonian H can immediately be carried out. The 



1) only the condition (43) plays a role. 



tilde on the integrated Hamilton function (41), however, is intended to emphasize that, 



even if the equations of motion after the normal form expansion can be derived from it 
according to Hamilton's equations, it will, in general, not represent the energy of the 
system. 



In order to achieve both, i.e. the satisfaction of the conditions of integrability (44) 



(45 ) for systems with d> 2 degrees of freedom and the local equivalence of the integrated 



Hamiltonian with the energy functional, a final transformation q q and p ^ p 
is necessary. The latter makes use of the fact that the transformations performed 
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in sections 2.3.1 and 2.3.2 still leave some freedom, e.g. the choice of multiples of 
eigenvectors when diagonalizing the equations of motion. 

For this purpose, we scale the generalized coordinates q ^ q and the momenta 
p — )■ p in (35)-(36) and (40) with some time-mdependent functions i>q-{q,p) and 



{q, p) according to 



Qi = J^q, {q, P) Qi , Pi = ^p, {q, P) Pi- (46) 
Since the qi and pi do only occur as products piQi in the integrated Hamiltonian H as 
well as in the transformed energy functional (40), the precise form of the Ug., Vp. is not 
of interest. We only demand their product 



(47) 



also to be a power series of the products qjPj which is necessary in order to preserve the 
structure of the respective equations. 

Because both sides of the equations (35) and (36) are scaled according to equation 
(46), one of the functions Vq^q^p) and Up-{q,p), respectively, cancels out, and the 
remaining products Vq^Vp^ can be replaced by the expansion (47). After rearranging 
the right-hand sides of the respective equations the coefficients of the expansions then 
depend on the and the conditions of integrability (44)-(45) provide a linear system 
of equations for their determination. The latter, however, is underdetermined and to 
obtain a complete set of equations one also requires the equivalence of the integrated 
Hamiltonian and the transformed energy functional. 

The cancelling out of one VqXq^p) or ^pAq^p) in the equations of motion leads to 
the crucial result that the occurrence of the fim in the integrated Hamiltonian H and 
in the energy functional exactly differs by one order, 

H= Ho + Hi{q,p) + H2{q,P, IJm,\m\=2) + H3{q,p, IJ,m,\m\<3) + ■ ■ ■ , (48) 
E= + El(qr,p,yU^,|m|=2) + ^2(g,P,/Um,|m|<3) + ^3(g,P,yUm,|m|<4) + • • • , (49) 

i.e. the order n term En of the transformed energy functional contains those fim with 
\m\ < n + 1, whereas in the corresponding term of the integrated Hamiltonian H only 
the scaling coefficients firn with \m\ < n occur. 

Consequently, orderwise comparing the coefficients of the expansions (48)-(49), 
yields both the constant of integration Hq and the required equations to uniquely 
determine the fim- The latter are, finally, determined orderwise as the solutions of a 
linear system of equations considering that the conditions of integrability (44)-(45) and 
the requirement of the equations (48)-(49) to be identical are fulfilled simultaneously. 

After this proper choice of the /im and the definition of action-variables by Jj = qiPi 
if Ji corresponds to a real eigenvalue Aj and Jj = iqiPi in case of purely imaginary ones 
icuj, we have 

H (J) = E{J) = H (J) . (50) 

This serves as classical Hamilton function in the sense that it reproduces the energy of 
the system, given by equation (15), and at the same time its Hamiltonian equations of 
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motion equivalently to equation ( 12 ) describe the dynamics in the vicinity of the fixed 



point up to some previously chosen order. 



2.4- Thermal decay rates 



After the classical Hamilton function (50) has been constructed, classical TST can be 
applied. In this section, we review the formulas to calculate the thermal decay rates 
which follow from this theory. 

We now assume the Hamilton function in phase space to exhibit a local minimum 
corresponding to its metastable ground state. Moreover let there exist an equilibrium 
point of saddle-centre-. . .-centre type which is the only channel the system may decay. 

The decay rate at a fixed energy is then given by the fiux through a dividing surface 
which divides the phase space into a region of "reactants" and "products" , respectively 
(see reference Il2]). If the system is in contact to a bath of finite temperature, the 



thermal decay rate is the Boltzmann average of this fiux. 

2.4- 1- The classical case The directional fiux through the dividing surface between the 



reactant's and product's region in phase space for a fixed energy is 12, 17, 18 



fiE) = {2'nY'^V{E), (51) 

with V{E) being the phase space volume of the actions {J2, ■ ■ ■ , Jd) which is enclosed 
by the contour ifcNF(0, J2, ■ ■ ■ , Jd) < E and Ji = corresponding to the "unstable 
direction" of the saddle. The thermal decay rate is then given by the Boltzmann average 



of equation (51) which, after a short calculation, yields 



Tci = o / e-''^cNP(o,J.,...,J,)dj^ _ _ _ ^j^^ (52) 



where (3 = l/k-BT, and Zq is the canonical partition function (cf. reference |19|). 
Since nearly all states will be localized in the vicinity of the local minimum, we will 
approximate the latter by 

with Hq^p{J[, . . . , J^) being the normal form expansion at the local minimum. 



dJi . . . dJd e-^^cNF(-^iv-^^) (53) 



In general, both integrals in equations (52) and (53) will not converge, which is due 



to the fact that, when actually computing the normal form, the expansion has to be 
truncated at some order. To compensate this, we restrict the phase space volume over 
which is integrated to the condition 

uji = —Q-j— > 0, z = 1, . . . , d (54) 

taking into account the fact that all frequencies occurring on the tori in phase space 
have to be non-negative. 
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2.4-2. The quantum case For the quantum mechanical case in equation (|8j), we proceed 
analogously. The quantum mechanical reaction probability for a one-degree-of-freedom 



system at a fixed energy is given by 12 , 14 



where the action variable J is implicitly given by Hq^p{J) = E. 

Again the thermal decay rate Fqm is given by the Boltzmann average of equation 



(55) which, after integration by parts, reads 

1 r g-/3J/QNF (J)-27r j(£;)/ri 
^^"^ " W(3Z~o J (i+e-2-J{i?)A)2 • (56) 
In the quantum case the canonical partition function is evaluated as a sum 

Zo = ^e-^^W(^'='-*("+V2))^ (57) 

n 

and because of the same reasons as in the classical case, we restrict the phase space 



volume for the integration in equation (56) and the summation in equation (57) to the 



condition (54) with i^cNF — > -f^QNF- 



3. Application and results 

3.1. The cubic model potential 

To demonstrate the applicability of the procedure to calculating thermal decay rates 
within the framework of a variational approach to the quantum mechanical wave 
function, we consider a cubic model potential 

V(x) = -^x^(2x-3-f). (58) 

It features a local minimum at Min(0, 0) and a local maximum at Max(7, a) for a, 7 > 
(see figure [T]) . 

Within the variational approach, we approximate the particle's wave function by a 
frozen Gaussian 

ip{x,t) = Afexp (— a(s — + ip{x — g)) (59) 
= Afexp (— ax^ + zx) (60) 

with fixed width a > 0. The centre q and momentum p of the wave function are 
combined in the complex variational parameter z with real and imaginary part Zj. = 2aq 
and Zi = p, respectively, and A/" = (7r/2a)"'^^^ exp (z'^/Aa) = A/'exp {—aq^ — ipq) gives the 
normalization of the wave function: J dx |^/'(a;,t)|^ = 1. With this definition every point 
in the complex z-plane can be identified with a certain position and state of motion of 
the Gaussian wave packet and vice versa. 

In this paper, we use a single frozen Gaussian, since in this case the energy functional 



(15) and the equations of motion for the variational parameters (12) can be obtained 



explicitly. Note, however, that the procedure is neither limited to Gaussian wave 
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Figure 1. Schematic drawing of the cubic model potential for a,j > (dashed line). 
The local minimum at Min(0, 0) is the metastable ground state of the system. A 
particle in the ground state can, after sufficient thermal excitation, cross the barrier 
formed by the local maximum at Max(7, a) and will then escape to x — oo. Within 
the variational approach, the particle is described by a Gaussian wave function of fixed 
width (solid line), and its two stationary states are located at the minimum and the 
saddle. 



functions nor to the use of a single Gaussian, but can without any change also be 
apphed to wave packets of coupled Gaussians (cf. reference [o]), or of different form. 

The energy functional of the system is given by the expectation value of the 
Hamiltonian [h = m = 1) 



H 



with the potential V{x) in equation (58). For the wave function (60) it reads 



E 



a (3a^7 + 30^^.(7^^ — 1) 



(61) 



(62) 



The application of the TDVP gives the equations of motion for the time evolution of 
the variational parameter and, after splitting into real and imaginary part, yields 

(63) 
(64) 



?>azl 



3aZr 



3a 



20273 



07^ 



2073 



TST can be applied to calculating the reaction rate of this system because of the 
following reason: First, we define "reactants" as initial states of motion for which the 
wave packet remains in the potential well, i.e. x < 7 for t — 00, and "products" as 
those which escape to a; — )■ cxo for t 00. With this definition the complex z-plane can 



be divided, by solving the equations of motion (63)-(64), in a region of reactants and 



products, respectively. The reaction rate is then given by the flux through the dividing 
surface between these two regions. 

Thermal excitation now enables a wave packet corresponding to the reactants to 
cross the barrier of the potential V{x) and to become a product. For low temperatures, 
the reaction path will cross the dividing surface at its point with the least energy. This 



point in the ^-plane is a saddle of the energy functional (62), consequently the decay 
rate of the system is given by the flux over this saddle. 
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However, since z is a "bad" variable in the sense that - as aheady stated in section 



2.3 - the equations of motion (63)-(64) cannot be derived from the energy functional 
by the canonical equations, the latter cannot be regarded as a Hamilton function which 
allows for the calculation of the flux. 

In order to obtain a Hamilton function H{J) which can be used for this purpose, 
we perform the steps described in section |2.3[ i.e. we determine the fixed points of the 



equations (63)-(64), expand them in their vicinity and order by order perform a normal 



form transformation. After integrating and adapting the result to the transformed 



energy functional, we end up with the desired Hamiltonian (50) which is locally 



equivalent to the equations (62)-(64). To illustrate the procedure in a more descriptive 
way, we provide a numerical example in the appendix for a certain set of parameters. 
Note that - as it is also the case in the analogous classical system of a point particle 



in the potential V{x) - the equations of motion (63)-(64) exhibit two fixed points. One 
of them corresponds to the metastable ground state of the system and the other to an 
unstable excited state, whose corresponding wave function is located in the vicinity of the 
potential's local maximum. While the ground state also exists in the framework of the 
exact Schrodinger equation, the existence of the unstable excited state is a consequence 
of the ansatz of the wave function within the variational approach. A wave packet which 
is located at the maximum of the potential V{x) is, of course, not a stationary solution 
of the exact Schrodinger equation, but would dissolve. Within the variational approach, 
this dissolving character, however, is taken into account by the instability of this state. 

Note furthermore, that the eigenvalues of the linearised equations of motion at 
the fixed points, which crucially determine the normal form expansion and, with it, 
the constructed Hamiltonian H{J), approximately reproduce the eigenfrequency of the 
harmonic oscillator which one obtains by harmonically approximating the potential V{x) 
at its minimum and maximum, respectively. For wave functions with smaller extension 
(larger a) the eigenvalues converge to the oscillation frequency of the harmonic oscillator 
and in the limit a — )• oo exactly reproduce it as is expected, since the dynamics of a 
strongly localized wave packet will be the same as that of a point particle. 

3.2. Comparison of the decay rates 

In order to compare the thermal decay rates calculated by the different methods we 
first determine the quantum and the classical normal form of a point particle in the 
model potential V{x) from equation ([s]) and its normal form analogue of the variational 



parameter space as described in section 2.3 The thermal decay rates are then given 



by equation (52) for the classical case and equation (56) for the quantum case (for the 
latter we set h = 1 throughout). 

Figure |2] shows a comparison of the thermal decay rates calculated for the cubic 
potential by the three methods in dependence of the parameter 7 of the potential and 
for different width parameters a of the Gaussian trial wave function. The temperature 



is set to /3 = 1.5 and we use a barrier height of a = 10 in equation (58). 
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Figure 2. Comparison of tlie tliermal decay rates of a particle placed in the metastable 
ground state of the cubic potential calculated by the classical normal form (solid line) , 
the quantum normal form (dashed line) and the variational approach with a frozen 
Gaussian wave function of different width a (dots). The data have been calculated 
in dependence of the parameter 7, for a temperature of /3 = l/k^T ~ 1.5 and a 
barrier height of a = 10 in equation (58 1. In the limit a — >■ 00 of a very narrow wave 



function (empty circles) the variational ansatz covers the results of the classical normal 
form while the decay rates obtained from this method increase with decreasing width 
parameter a (triangles and filled circles). For a broad wave function (squares) it shows 
a perfect match with the decay rate obtained from the quantum normal form of the 
corresponding point particle. 



The classical (solid line) and quantum (dashed line) normal forms reveal a 
significantly differing decay rate which is especially pronounced for small 7 due to the 
quantization effects in the narrow potential well. The variational approach with a single 
Gaussian trial wave function, however, is able to reproduce both curves. In the limit 
of a very narrow wave function (empty circles), i.e. a — )■ 00, we find a perfect match 
of the latter with the classical result. When we increase the width of the Gaussian 
wave function, also the decay rate increases (triangles and filled circles). For a width of 
a = 1.9 (squares), we finally recover the decay rate calculated from the quantum normal 
form of the corresponding point particle. 

The classical limit of the decay rate can easily be understood considering the 
classical limit of the variational wave function. For the narrow wave function, only 
the local properties of the potential will be of importance so that the potential can be 
approximated harmonically in this limit. The characteristic width of the Gaussian is 
then obtained from that of the ground state wave function of the harmonic oscillator, 
which is given by a = u/{2h), with u = ^/\V"{x)\^^q. The classical limit h ^ then 
yields a — )■ 00 for which we observe convergence of the decay rate obtained from the 
variational approach to the one obtained from the classical normal form. For broad 
wave functions the anharmonicity of the potential becomes relevant and the width of 
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Figure 3. Thermal decay rates for the cubic potential for different temperatures. 
As in figure [2] we compare the results of the classical (solid lines) and quantum 
(dashed lines) normal forms as well as the variational approach in the limit of broad 
(squares) and narrow (empty circles) Gaussian trial wave functions. With increasing 
temperature, also the decay rates increase and the difference between the classical and 
quantum normal forms becomes smaller in the same way as the two cases of narrow 
and broad wave function do within the variational approach. For /3 < 0.1 the decay 
rates calculated by the different methods cannot be distinguished any more. 



the harmonic oscillator will be only a crude approximation. For h = 1 the width of 
a = 1.9, for which we observed good agreement with the quantum normal form result, 
is in accordance with the width of the harmonic oscillator ground state within about a 
factor of 2 over the whole range. According to this, intermediate values 1.9 < a < oo 
correspond to effective values < < 1 of the Planck constant. 

We further observe this behaviour at different temperatures as can be seen in figure 
|3] for some exemplary values. Again, the results using the classical normal form (solid 
lines) are perfectly recovered within the limit a — i- oo of very narrow wave functions 
(empty circles) and for broad ones (squares) the decay rates agree with those obtained 
from quantum normal form (dashed lines). 

With increasing temperature, the calculated decay rates also increase rapidly by 
several orders of magnitude. However, the difference between the classical and the 
quantum normal forms becomes smaller as well as it does for the two limits of the wave 
function concerning their width. At high temperatures compared to the height of the 
potential barrier we reach the classical limit. For the parameters used here, this is the 



TST for wave packet dynamics. I. Thermal decay in metastable Schrddinger systemslS 



case for /3 < 0.1, where the difference of the decay rates calculated by the three methods 
vanishes and they cannot be distinguished any more. 

We emphasize that - in contrast to the classical and the quantum normal form - 
the method with wave packets described by the TDVP does not require a Hamilton 
function for the system to be known at the beginning of the procedure. Instead, only 
an expression for the energy functional and a set of differential equations describing the 
time evolution of the variational parameters are necessary. The required Hamiltonian 
is, after the mapping to phase space coordinates, a natural result of the transformations 
discussed in this paper. 



4. Summary and outlook 

We have presented a method to apply transition state theory to wave packet dynamics 
in metastable Schrodinger systems, which is based on local expansions of the energy 
functional and the equations of motions derived from a time-dependent variational 
principle in the vicinity of the barrier fixed point of saddle-centre-. . . -centre type. 
Fulfilling the integrability conditions, the equations of motion can, after a normal form 
expansion, be integrated to a common Hamilton function expressed in action-variables. 
Its equivalence to the energy functional is guaranteed by an additional transformation 
adapting it to the latter. 

To demonstrate the method, we calculated thermal decay rates of a particle placed 
at the metastable ground state of a cubic potential. On the one hand the decay rate was 
obtained by the Hamilton function of a point particle using the classical and quantum 
normal forms. On the other hand this was done by the presented method with a frozen 
Gaussian wave function of different width described in the framework of a variational 
approach to the Schrodinger equation and the results show an excellent agreement with 
the well established methods. 

The application of the variational approach to wave packets within the framework 
of a time- dependent variational principle together with the demonstrated mapping to 
classical phase space further allows to calculate decay rates for quantum mechanical 
systems where an analogous classical Hamiltonian is not directly accessible. 

Such systems are, e.g. Bose-Einstein condensates with additional long-range 
interaction which have been successfully approached variationally with coupled Gaussian 



wave functions 10 ,11 . In case of an attractive short-range contact interaction, their 
ground state is metastable and a sufficient thermal excitation may lead to the collapse 
of the condensate by crossing a barrier which is formed by an unstable excited state 
that reveals the properties of a saddle-centre-. . .-centre equilibrium point. Thus, the 
variational approach to BECs is predestined for the procedure presented in this paper 
and we refer the reader to reference |9] for the application to these ultra-cold gases. 

Of course, this method can also be applied to the field of chemical reactions 
where the behaviour of reactants can be studied without approximations such as in 
the framework of an ab initio calculated Born-Oppenheimer potential energy surface if 
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a suitable parametrization of the system's wave function has been found. 
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Appendix: Example 



In order to give a descriptive illustration of the variational parameters' mapping to 
classical phase space, we present in this appendix a numerical example for the parameters 
a = 1.9, a = 10 and 7 = 2 in fifth order of the equations of motion and sixth order of 
the energy functional, respectively. 



For this set of parameters the Taylor expansion (19) of the equations of motion 



(63)-(64) in the vicinity of the stable fixed point reads 
Xi = 3.8x2, 

X2 = -3.67852x1 + 0.519391x2 



and for the corresponding expansion (37) of the energy functional (62) one obtains 
E = 1.90363 + 0.484015X? - 0.0455606x? + 0.5x^. 



(65) 
(66) 

(67) 



At this point, there is no obvious way to derive equations (65) and (66) from equation 



(67) via Hamilton's equations. 

With the eigenvalues iwi 2 = ±3.738771, the normal form expansion of the equations 
of motion and the introduction of the canonical variables q = Xi, p = X2 results in 

(68) 
(69) 



Q 
P 



3.73877ig - 0.063123ipg^ - 0.00250446ip2g^ 
-3.73877ip + 0.063123ip2g + 0.00250446ip^g2 



for equations (35) and (36), where the higher-order terms have arisen during the 



transformations. The coefficients on the right-hand side of these two equations only 



differ in sign and therefore fulfil the conditions of integrability (43)-(45). Consequently, 
they can easily be integrated to a Hamilton function 

H = Ho + 3.73877ipg - 0.0315615i(M)^ - 0.0008482i(pg)^ (70) 

However, after having applied the change of coordinates corresponding to the 



transformation above, the energy functional (40) reads 

E = 1.90363 + 0.983756pg - 0.00830456(pg)2 



0.000219661(pg) (71) 



whose coefficients do not agree with equation (70). 



To achieve this, we apply the scaling of the phase space variables q ^ q and p ^ p 



according to equations (46)-(47) which after integration of the equations of motion 



TST for wave packet dynamics. I. Thermal decay in metastable Schrddinger systems2Q 



together with the definition of the action-variable J = ipq leads to 

H= #o- 3.73877J + 0.0315615iAiiJ^ 

- (0.000834821/i^ + 0.021041/X2) J^ (72) 



E = 1.90363 + 0.983756i/ii J + (0.00830456^^2 - 0.983756^^2) 

+ (0.000219661i/i? + 0.0166091i/ii/i2) (73) 

Setting the constant of integration to i^o = 1-90363 and choosing /ii = 3. 80051, /i2 = 0.0 
finally yields 

H{J) = H{J) = E{J) 

= 1.90363 - 3.73877J- 0.119949J2 + 0.012058J^ (74) 

which is the desired Hamilton function in third order of the action-variables. 

We note that the vanishing scaling parameter fi2 = 0.0 is a consequence of the 
ansatz (60) where the real and imaginary part of the variational parameter z already 
have the meaning of position and momentum of the wave function. In general, all the 
firn will be non-zero. 
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